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Mesoscopic molecular dynamics simulations are used to determine the large scale structure of 
several binary polymer mixtures of various chemical architecture, concentration, and thermody- 
namic conditions. By implementing an analytical formalism, which is based on the solution to the 
Ornstein-Zernike equation, each polymer chain is mapped onto the level of a single soft colloid. 
From the appropriate closure relation, the effective, soft-core potential between coarse-grained 
units is obtained and used as input to our mesoscale simulations. The potential derived in this 
manner is analytical and explicitly parameter dependent, making it general and transferable to 
numerous systems of interest. From computer simulations performed under various thermodynamic 
conditions the structure of the polymer mixture, through pair correlation functions, is determined 
over the entire miscible region of the phase diagram. In the athermal regime mesoscale simulations 
exhibit quantitative agreement with united atom simulations. Furthermore, they also provide 
information at larger scales than can be attained by united atom simulations and in the thermal 
regime approaching the phase transition. 
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I. INTRODUCTION 

The mixing of two or more types of polymers is of 
great scientific and technological interest, as it opens 
up the possibility of creating new materials emerging 
with specific physical and chemical properties. |T] How- 
ever, although polymer blends have been very much a 
part of everyday life, these continue to be a source of 
extensive scientific inquiry. The rich physics in polymer 
mixtures stems in part from the fact that their struc- 
ture and dynamics change as thermodynamic conditions 
that cause phase separation (i.e. spinodal decomposition) 
are approached. Mixture stability is not only driven by 
temperature and composition, but also by differences in 
chain length and monomer architecture that may con- 
tribute substantial entropic effects. From physical and 
engineering standpoints, a goal is to understand and pre- 
dict changes that a polymer will undergo when mixed 
with another polymer system [21 [3] . Theoretical efforts 
have focused on this area, resulting in the development 
of several approaches connecting monomeric structure to 
static and dynamic properties as a function of thermo- 
dynamic parameters |3HS]- 

Molecular dynamics (MD) studies have aided our un- 
derstanding of the correlation between local (intra- and 
intermolccular) structure and global fluid properties that 
govern polymer mixtures |10H16| . However, these compu- 
tational approaches have been limited to relatively small 
length and time scales because present-day computing 
power limits the attainment of extended scales. [17]. For 
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polymer blends, this is particularly problematic because 
as the spinodal decomposition is approached, a diver- 
gent length scale in concentration fluctuations develops, 
thereby readily exceeding box sizes commonly used to 
model polymeric ensembles. In this manner, the deter- 
mination of large-scale properties for polymer mixtures 
with atomic resolution requires increasingly larger simu- 
lation boxes, and becomes rapidly unfeasable. 

One way to circumvent this difficulty is to introduce a 
coarse-grained description of the polymer mixture, where 
small-scale atomistic details are statistically averaged 
out. The resulting structure allows for computer sim- 
ulations of large systems, including a high number of 
polymers, chains with large molecular weights, and for a 
long timescale. Recent mesoscopic descriptions based on 
various coarse-grained approaches have been presented 
in the literature and effectively utilized to investigate 
the demixing of polymer solutions |18j. immiscible poly- 
mer blends[19 , and star-linear polymer mixtures. |20| 
By combining information obtained from simulations of 
the coarse-grained system with information obtained at 
smaller lengthscales and shorter times, in the context of 
a multiscale procedure, it is possible to obtain a com- 
plete and exhaustive description of the polymer mixture 
at all lengthscales of interest. [21] We have previously 
shown that a procedure, which combines information at 
the two length scales of the effective segment length and 
the radius of gyration, accurately determines the struc- 
tural properties of homopolymer melts. [22] This multi- 
scale modeling procedure has been recently extended to 
treat polymer mixtures. [23| 

In this publication, we focus on performing simulations 
of polymer mixtures where the ensemble of polymers is 
mapped onto a system of interacting soft, colloidal par- 
ticles. At this level of coarse-graining each polymer is 
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represented as one soft colloidal particle centered at the 
polymer center-of-mass. The concept of modeling flexi- 
ble polymer chains as soft spheres was originally proposed 
by Flory and KrigbaumfM| for dilute polymer solutions; 
however, the resulting effective potential exhibited an 
incorrect trend in repulsive interactions with increasing 
chain length or polymer density. For dilute to semidi- 
lute solutions, Louis and co-workers were able to obtain 
correct scaling behavior for polymer chains mapped onto 
soft spheres using liquid state theory in conjunction with 
Monte Carlo simulations. [5S] This approach was later ex- 
tended to poor solvent conditions where coils contract to 
avoid contact with the solvent. Importantly, this work 
shows that for a given set of state conditions, a simple 
effective pair potential, v{r), is capable of reproducing 
the radial distribution function, g{r), which depends on 
many body interactions. Due to this reduction in inter- 
actions that need to be calculated at each computational 
step, mesoscale simulations of polymers represented as 
soft colloids, are useful in determining many structural 
as well as dynamical properties of polymer blends, includ- 
ing how morphologies develop depending on the thermo- 
dynamic parameters of the system [27]. The advantage 
of this extreme coarse-graining of the polymer is that 
it is possible to simulate very large ensembles of parti- 
cles, i.e. adopt large simulation boxes, with a modest in- 
crease in the computational time. Such an extreme level 
of coarse-graining becomes important for simulations of 
systems where the relevant range of length scales is par- 
ticularly large, for example in micellar aggregates of ionic 
surfactants. [55] 

The implementation of mesoscale simulations requires 
attaining the effective "bare" potential, w'^'^(r), that char- 
acterizes pairwise decomposable interactions between 
molecular center-of-mass (com) sites. Since w^^(r) cor- 
responds to a free energy obtained from the monomer 
frame of reference, it depends on the thermodynamic 
state of the system, as specified by the density, tempera- 
ture, polymer molecular weight, and composition of the 
mixture. 

The determination of a reliable, fully transferable, 
coarse-grained potential is a highly desirable goal.[^ 
While phenomenological as well as rigorously numerical 
approaches have been described in the literature [25l l30l - 
132] to determine effective coarse-grained pair potentials 
for polymer systems, their reliance on acquired micro- 
scopic simulation data partially defeats the gain in com- 
putational time that is possible with a coarse-graining 
procedure. Typically, mesoscopic potentials are opti- 
mized to full atomistic simuluations under a given set of 
thermodynamic conditions, limiting their transferability. 
Moreover, the optimization procedure generally starts 
from the mean-force potential, which is the effective po- 
tential between two particles in the field of the surround- 
ing particles, and is conceptually different from the bare 
potential, w'^^(r). We have developed an analytical for- 
malism to calculate the effective potential, v'^'^{r), start- 
ing from liquid state theory and solving the Ornstein- 



Zernike equation. [53] [M] The potential obtained in such 
a manner is explicitly related to the structural parame- 
ters of the polymer, making it specific to any system we 
desire to simulate, but also fully transferable to systems 
with different molecular structure and thermodynamic 
conditions. 

More specifically, our non-phenomenological expres- 
sions for the com-com total pair intermolecular correla- 
tion functions, ft-^'^^(f), between self (aa) and cross (a/?) 
interactions, are obtained from a generalized Ornstein- 
Zernike integral equation for binary polymer mixtures 
where atomistic sites are accounted for as real sites, 
while coarse-grained sites are treated as auxiliary sites. 
The equation formally bridges information from the mi- 
croscopic (monomer) domain to mesoscopic (molecular) 
scales. The derived equations for h'^^^r) reproduce well 
and with no fitting parameters united atom (UA) molec- 
ular dynamics (MD) simulation data in both real and 
reciprocal spaces. Our renormalized description, [321 [31] 
correctly recovers the known analytical expressions for 
density and concentration fluctuations of mixtures of col- 
loidal liquids, by Kirkwood and Buff, [35^ and Bhatia and 
Thornton, [32] slightly modified as our expressions ac- 
count for soft potentials instead of hard-core ones. To- 
gether, these tests provide a benchmark which supports 
the foundation of our renormalization procedure. 

In the current publication we extend our previous work 
to demonstrate that the derived coarse-grained effective 
potential, when input to mesoscale MD simulations of bi- 
nary polymer blends, is capable of reproducing the large 
scale structure of the polymer, which depends on many 
body interactions at the monomer level. We demonstrate 
that our approach is useful to produce mesoscale simu- 
lations of binary homopolymer mixtures at various tem- 
peratures and concentrations approaching the spinodal. 
In this way, we are able to test the derived potential in 
terms of how well it reproduces the mesoscopic struc- 
ture of the mixture over different regions of the phase 
diagram, where atomistic level simulations are compu- 
tationally exhaustive. We make use of the hypernetted- 
chain (HNC) closure to calculate the effective potential, 
v'^l^{r), as a function of the total correlation function, 
h'^p{r). In turn, the pair potentials, V^i^ir), are used as 
an input to a coarse-grained simulation, where polymers 
interact as soft colloidal particles. 

Systems investigated are blends of polyethylene, poly- 
isobutylene, as well as polypropylenes in their head-to- 
head, istotactic, as well as syndiotactic forms. We show 
that our method is robust and allows for the equilib- 
rium structural properties of the fluid mixture to be 
readily calculated under any thermodynamic parameters 
of interest. While the focus here is on static proper- 
ties, the derived potential is widely applicable to non- 
equillibrium systems, and may be useful in other meth- 
ods commonly employed, such as dissipative particle 
dynamics [23, which currently rely on numerical poten- 
tials. 

To demonstrate the reliability of our mesoscale sim- 
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ulations, we first consider systems which are far from 
the spinodal temperature, for which the Uquid is well- 
mixed, and show that in this regime, our ^^'^(t') repro- 
duces quantitatively the liquid structure at the level of 
com obtained from UA MD simulat ions 1 1 3 1 FH] . establish- 
ing consistency between the different levels of description. 
We then extend our analysis to different thermodynamic 
conditions for which no UA data is available, focusing 
on mixtures of head-to- head polypropylene and polyethy- 
lene (hhPP/PE) in the miscible region of the phase dia- 
gram. The temperature dependence is expressed in terms 
of a single interaction parameter that enters our simula- 
tion through the analytical form of the potential and de- 
pends on the local interactions between monomers. From 
the simulation trajectories, we calculate pair distribution 
functions which depend on the parameters of the sys- 
tem and manifest the trends for demixing of the coarse- 
grained mixture. From this treatment, and from the fact 
that our expressions recover known expressions for col- 
loidal liquids, we calculate the concentration fluctuation 
contribution to the scattering function, S{k), which is 
related to the scattered light intensity measured in ex- 
perimental studies of critical binary polymer mixtures. 
In this way, we are able to test the derived potential in 
terms of how well it reproduces the mesoscopic struc- 
ture of the mixture over different regions of the phase 
diagram. Results show that the effective potential be- 
tween center of masses of polymers in a mixture can be 
used to produce mesoscopic simulations under a variety 
of thermodynamic conditions, making the procedure use- 
ful to a number of application in polymer physics. To 
test the versatility of our approach we show how the pro- 
posed theory may be applied to blends that present a 
lower critical solution temperature (LCST). We study the 
hhPP /PIB mixture using the temperature dependence of 
the X parameter from the literature and run mesoscale 
simulations at various temperatures. The concentration 
fluctuation structure factor shows that fluctuations in 
concentration become smaller as the temperatures de- 
creases, as it is expected. 

The current publication is organized as follows. In 
Section Kll a review of our derivation for h^p{r) is pro- 
vided. Tnese results are then used in Section IIIII for 



the calculation of v^p{r), obtained from the HNC clo- 
sure. In Section 
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our representations for ^^'^(t') are 
implemented in mesoscopic simulations of binary poly- 
mer blends, where the mixture is modeled as an assem- 
bly of soft, interacting colloidal particles. Section |V] com- 
pares predicted total correlation functions, analytical and 
from mesoscale simulations, with data from united atom 



molecular dynamic simulations. In Section VI data ob- 
tained from mesoscopic simulations are used to calculated 
liquid partial structure factors which are used to deter- 
mine the phase diagram of the mixture. Section |VlT] in- 
vestigates how corrections to the Debye form factor affect 
the precision of the predicted total correlation functions. 
Finally, in Section |VlTl] we show how our approach is ef- 
fective in predicting high temperature demixing for the 



LCST blend hhPP/PIB. 



II. MESOSCOPIC PAIR CORRELATION 
FUNCTIONS FOR ASYMMETRIC BINARY 
POLYMER BLENDS 

Pair correlation functions (pcfs) comprise a standard 
approach in liquid-state theory to describe the structure 
of liquids. Since it is generally sufficient to account for 
two-body correlations, pcfs can be employed to deter- 
mine all structural and thermodynamic properties of a 
system [38 . Moreover, these are input to the equation of 
motion describing cooperative dynamics, where the effect 
of the surrounding medium on single-chain dynamics is 
taken into account ^39^ ,40, . In the context of our coarse- 
graining procedure, pcfs provide a convenient manner of 
calculating, with the aid of an appropriate closure, the 
bare pair potential needed to perform mesoscopic simu- 
lations. In this section, a brief review is given of the the- 
oretical formalism we previously developed to describe 
coarse-grained binary mixtures of asymmetric polymer 
blends [33IMI- 

Our model for a binary blend consists of A and B ho- 
mopolymers, having Na and Nb monomer sites with seg- 
ment lengths a A and as, respectively. For simplicity, 
these monomer sites are taken to span equivalent vol- 
umes so that the polymer volume fraction is given by 
(j) = nANA/{nANA + ubNb), where Ua is the number 
of molecules of type a in the mixture with a G {A, B}. 
While p = {uaNa + nBNB)/V quantifies the total num- 
ber of monomer sites contained in a region of space 
spanned by V ^ the site and chain number densities for 
molecules of type A are given by pA — tiaNa/V — (f>p 
and pc,A = nA/V, respectively. 

The derivation of our analytical expressions for total 
intcrmolecular com pcfs in a polymer mixture extends 
from a procedure outlined by Krakoviack, et al. |41j for 
homopolymer solutions. The key step in this approach 
is to include molecular coms as auxiliary sites, along 
with monomer sites serving as real sites. The generalized 
Ornstein-Zernike integral equation is solved in reciprocal 
space and is given by 



H(/c) = n{k)c{k) [n{k) + H(fc)] 



(1) 



where H(fc) is the matrix of total intcrmolecular pcfs, 
C(fc) is the matrix of direct intcrmolecular pcfs, and n{k) 
represents the matrix of intramolecular pcfs. Specializ- 
ing to the case of a binary polymer mixture, each matrix 
in Equation [T] is of rank four, composed of four 2x2 
blocks that account for monomer- monomer (mm), com- 
com (cc), and the corresponding cross (cm and mc) in- 
teractions. For instance, 



H(fc) = 



(fc) 
(fc) 



H 



(fc) 
(fc) 



(2) 



The remaining matrices in Equation [llfollow an arrange- 
ment analogous to that of Equation |2j Each block in 
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Equation [2] contains self {aa) interactions along its diag- 
onal, whereas cross (a/3) interactions occupy off-diagonal 
positions. 

As a next step, the individual block elements 
that define the matrices in Equation [l] are defined. 
The intermolecular total pcf matrix H(fc) is com- 
posed of the chain-averaged monomer-monomer pcfs 
{k) = paPph^p^{k), the com-monomer pcfs 
[k) = pcaPph^^'pik) = HJ^^{k), and com-com 
pcf H^lik) = PcaPci^hafsik). Note that in gen- 
eral, h'^^ik) = hf^{k) while /i^'^(fc) ^ h^^ik) when 
a 7^ /?. The intramolecular pcf matrix Cl(k) is sim- 
ilarly composed of = PaP/3W™™(fc)(5a/3Za/3, 

^a^(fc) = Pc,c.Ppoj^'^{k)So.pz^p = n^^{k), and = 
PcaNpPc^pdapZap, wherc Zap = [(f>i3P{'2 - Safj)]^'^ ■ Fi- 
nally, the intermolecular direct total pcf matrix C{k) 
acquires a much simpler structure when including the 
approximation that any interaction involving auxiliary 
sites is negligible. Under these conditions, the only 
non-zero block in C(fc) involves monomer- monomer pcfs. 



defined as C^^^^lfc) 



Using the matrix definitions described above. Equation 
[TJis solved to obtain the intermolecular mesoscopic total 
pcfs, which are given by the relation 



uj^^{k)u:^^{k) 



(k). 



(3) 



Upon inspection, it is readily seen that Equation [3] for- 
mally connects com distribution functions to monomer- 
monomer intra- and intermolecular distribution func- 
tions. In this manner, one calculates mesoscale prop- 
erties from information on the local polymer scale. As 
mentioned before, this feature is relevant because prop- 
erties on the mesoscale ultimately depend on small-scale 
interactions. 

To obtain analytical solutions for h^p{k), a brief de- 
scription is given for each of the correlation functions 
entering into Equation |3] The com-monomer intramolec- 
ular pcf can be approximated in reciprocal space with a 
Gaussian distribution as 



ujZik)=Nae 



(4) 



with the molecular radius of gyration defined as Rga — 
(A^/6)^/^(Tq. On the other hand, the monomer-monomer 
intramolecular pcf is given by the Debye formula. 



2Na 



- 1 + enl 



ga 



k*Rga 



(5) 



For analytical convenience, however, it is costumary to 
approximate Equation[5]with its Fade approxiniant given 

bygg 



<:rik) 



Nr. 



1 + k^mj2 



(6) 



Although approximated, inclusion of Equation |6] allows 
for a convenient analytic expression for h'^'^{r) given be- 
low by Equation [9] which has been shown to give good 
agreement with simulations for the total pair correlation 
function in both real and reciprocal spaces ^SSj. We dis- 
cuss the implication for this approximation in detail be- 
low. In the current publication, we use both Equation [5] 
and Equation [6] for a;"""(fc), and compare the resulting 
mesoscopic h'^^k) from Equation |3] 

The respective monomer-monomer intermolecular to- 
tal pcfs used are taken from the thread limit of the Foly- 
mer Reference Interaction Site Model gl [S] (FRISM). 
The initial analytical treatment in the context of PRISM 
for polymer mixtures [44j has been extended by Yatsenko 
et al. [M' to account for chain asymmetry effects in the 
system. In this approach, a new parameter enters the for- 
malism, 7 = aBl<^A-, which defines the monomer asym- 
metry. 

While the thread model for polymer chains coarsely 
describes the liquid structure on local scales, it accurately 
captures the onset of the "correlation hole" effect at a 
length scale of Rg. Given that the spatial dimension 
of interest in our description is i?g, the thread limit of 
FRISM is an adequate representation for the intended 
purpose of the present work. The solutions are given by 



/i™™(r) = — ^ 



1 '^Al 



where 



V24</)(l-0)x.(l-x/Xs) 



(8) 



is the length scale governing concentration fluctuations, 
which diverges at the spinodal temperature. Here, x is a 
single interaction parameter that depends on the specific 
nearest neighbor pair energies between two A A, AB, or 
BB monomers, and is given by x = cab — (^aa + eB-B)/2. 
In a mesoscopic treatment which averages out the spe- 
cific monomer interactions, x is an input parameter cor- 
responding to the temperature dependance of a specific 
polymer architecture. From our definitions it clear that 
the quantity xl P is the analog of the Flory-Huggins inter- 
action parameter, and at the spinodal temperature x ~^ 
Xs, where Xs = [2iVA0]"' + [2iVs(l-<^)]"^ The quantity, 
(1— x/Xs) can be seen as a reduced temperature that indi- 
cates how far the system is from its spinodal temperature. 
Also in EquationlTl ^ca = -Rga/2^/^ is the length scale of 
the correlation hole while = Trptr^^/S -I- ^^^^ is the 
density correlation length scale with ct^o — (/)^ct^ -l-^ctcr?. 



(7) 
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This latter definition reintroduces finite-size effects, local 
semiflexibility, and branching that pertain to each com- 
ponent through a melt-like description. The effective seg- 
ment length scales are determined from the radius of gy- 
ration of each component polymer, through the relation 
Oc. = i6/N^y/^Rg. 

Inserting the definitions from Eqs. Q, Q, and ([T]) into 
([3]), the intermolecular total pcfs at the com level read 



'■BB 



(r) 
(r) 



1 - 
I" 

^AB 



ltBir)+l-'l'BB(.r) 



(9) 



hABir) = -lUr)+I'ABir), 



where I^fsi^) ^-nd lapi^) identify the concentration and 
density fluctuation contributions, respectively. We intro- 
duce here a compact notation with the function I^j^ir) 
defined as 



R„ 



X [e''/«^erfc ( # 
_e-''/?Aerfc ( ^ 



gag 



_ r73 



2R. 



and 



-P2 -PT- 



?c|3/3 



f2 



(10) 



(11) 



(12) 



where G iU^^p} and Cp 

gyration in the blend are defined such that 2R^ 



^aa + i??S=4CL^i,withec 



3/{npa\g). Radu of 

gap 



The development presented here is the required input 
to the derivation of the effective pair interaction poten- 
tials, a topic that will be addressed in the following sec- 
tion. 



III. THE EFFECTIVE SOFT CORE POTENTIAL 

Our theoretical approach is based on an integral equa- 
tion description of intermolecular pair correlation func- 
tions. A closure approximation is needed to connect 
these liquid structure functions to the effective pair in- 
teraction potentials which are required to perform, in 
our case, mesoscopic simulations of polymer mixtures 
mapped onto ensembles of soft colloidal particles. Since 
the fundamental units in our description interact through 



a soft-core potential, use is made of the hypernetted- 
chain (HNC) closure, which is known to be accurate for 
such systems. [IS] The relationship connecting the effec- 
tive pair interaction potential w^^(r) with pcfs is given 
by the HNC closure as 



{kBT)- 



/i-^(r)-ln [l + /i-^(r)] -C^(r), 

(13) 

where c'^p(r) is the direct pcf. Following the matrix defi- 
nitions presented in Section [ll] and taking our system to 
be a simple liquid comprised by soft colloidal particles, 
the direct pair correlation functions are defined by 



Cq^/3(^) 



1 

Pc,a 



{Pc,a + Pc,[i) |Scc(fc)| 



iPc,a + Pc,l3) \Scc{k)\ 



(14) 



where S*^^ and are the static structure factors 



for a binary mixture, and |Scc(fc)| = S'^j^{k)S'j^g{k) — 
[S'^g{k)] is the determinant of the mesoscopic static 
structure factor matrix. For a binary mixture these static 
structure factors are given by 



SAAik) - (t) + (l)^Pchh''^Aik) , 
SBBik) = l-q^+il-^)^Pchhj^Bik) 

SABik) = 0(l-0)pch/iTs(fc) , 



(15) 



where the total chain density, pch — p/N . By inserting 



Eqs. (p3|) and (|14|) into (|T3|), the v^^^ir) are obtained 



Since the potential is obtained from an inversion pro- 
cedure utilizing the HNC closure, the adequacy of this 
method is limited to systems for which the HNC is valid, 
which includes the weakly interacting soft particles mod- 
eled here. However, for systems of hard spheres or for 
low density ionic fluid models, such as the restricted 
primitive model (RPM), the HNC has been shown to 
be inaccurate [46], and more sophisticated closures are 
required. [32. ,47 

The analytical solution of our mesoscopic approach 
represents an advantage to previous work, where effective 
pair potentials are derived from simulations performed on 
a microscopic level. Such a requirement partially defeats 
the computational time gains behind a coarse-graining 
procedure. Overall, our approach bypasses the need to 
perform atomistic simulations for each thermodynamic 
state point of interest, which is necessary in numerical 
implementations since the effective pair interaction po- 
tentials depend on the state of the system. This can be 
readily appreciated from the pcfs that enter into the HNC 
closure, which are themselves state-dependent. 

We investigated the effect that the use of the Debye 
formalism, Equationjsj or of its Fade approximant, Equa- 
tion[6j for the monomer form factor in the denominator of 
Equation [3] has on the calculation of the potential. The 
Fade approximant is less precise than the Debye equa- 
tion, but it allows for the analytical solution of the total 
correlation functions, Equation [9j We observe that when 



6 



Equation [6] is used, singular points arise in the low k 
regime in the solution of Equation 14 for c'^'^(k) as the 



determinant of the mesoscopic static structure factor, 
|S,,(fc)| = S'XA{k)S'j^B{k) - passes through 

zero. This corresponds to an unphysical region of neg- 
ative compressibility. When Equation [5] is used instead, 
such singular points do not arise. 

For homopolymer melts, it has previously been de- 
termined that the singularities in c'^'^(fc) occur as a re- 
sult of the intrinsic error introduced in Equation |6] by 
the Fade approximation. In order to obtain a us- 
able form of the effective potential from Equation [TSj 
we tested two schemes: in scheme 1, we enforced the 
condition that d'^ik = 0) < c='=(fc) < for low k, 
which effectively eliminates any singularities from the di- 
rect correlation function; in scheme 2, we enforced the 
isothermal compressibility limit, such that for regions 
where |Scc(fc)| < |Scc(0)|, we truncated h'^'^{k) so that 
^Q/3(^) ~ ^'api^ ~ 0). The two schemes are equivalent 
and give identical results within the precision of our cal- 
culation. This is so, because polymer liquids are almost 
incompressible . 

In this work, we study polymer blends of polyethylene 
(FE), polyisobutylene (FIB), and polypropylenes in their 
head-to-head (hhFF), isotactic (IFF), and syndiotactic 

While the full Debye form (Equation [s]) for the 
monomer form factor prevents an explicit analytic ex- 
pression for h'^'^{r) in the form of Equation [9j which was 
the motivation for adopting the Fade approximation, a 
numerically obtained h'^'^{r) still can be readily obtained 
for any given system and so does not represent a limita- 
tion to our approach and avoids any singularities in the 
low k region for c'^'^(r). For this reason, in our following 
calculations the Debye approximation will be preferen- 
tially used. 

The potentials ^^^(f), calculated following the proce- 
dures discussed here, are required to carry out the simu- 
lations of the polymer liquid on a mesoscopic level. In the 
next section, we discuss the implementation of the w^"^ (r) 
to our mesoscopic simulations and in the following sec- 
tions we compare mesoscale simulation results with UA 
MD simulations and theoretical predictions. 



IV. MESOSCOPIC SIMULATIONS OF BINARY 
MIXTURES 

In this section, we implement molecular dynamics sim- 
ulations for the systems presented in Table |lj and we 
describe our methodology for carrying out mesoscopic 
simulations. 

In our mesoscopic modeling approach, we implement 
classical MD simulations, where the ensemble is evolved 
in the microcanonical [N, V, E) ensemble. In the ini- 
tialization step, all particles are placed on a cubic lattice 
with periodic boundary conditions. We use reduced units 
such that all length units are scaled by Rg {r* = r/Rg) 
and energies scaled by fc^T. The type of particle that 



(sFF) forms. The effective pair potential, u°'^(r) for in- 
teractions of type AA, BB and AB is calculated for the 
different binary polymer mixtures and for hhFF:FE un- 
der different values of (p and x using both the Debye form 
and Fade form of the intramolecular distribution function 
(Equation [5] and |6]). As a model calculation of the poten- 
tial, we present the results for the prototypical hhFF /FE 
polymer blend in Figure[T] which shows how the potential 
depends on the reduced temperature (1 — xlXs)- Input 
parameters to our theoretical calculations are reported in 
Table |l] as data for the U A simulations against which we 
test our approach. [T31 [Tl] Although there is a noticeable 
difference in the potential obtained using either Equation 
[5] or [6] they are qualitatively similar in many respects. 
For example, under athermal conditions, the mixture is 
random and the number of AB contacts is in between 
those of the self terms, AA and BB. Correspondingly, 
pair interactions accounting for AB contacts must be in- 
termediately repulsive. This effect is reflected in the plot 
of v^p{r). The A- type (flexible hhFF) particles display 
the highest repulsive response as a consequence of their 
stronger correlation hole effect. The inset of Figure [T] 
highlights the change in the repulsive component in the 
potential, as the ratio, x/Xs, is varied. 

TABLE I: Polyolefin blends (T = 453 K, Na ^ Ng = 96) 



Blend[A/B] 


p [sites/ A^] 


R,A [A] 


7 


X 


hhPP/PE 0.50 


0.0332 


12.32 


1.34 


-0.0294-1- 17.58/r"''' 


PIB/PE 0.50 


0.0343 


9.76 


1.68 


0.00257 4- 4.99/T' 


PIB/sPP 0.50 


0.0343 


9.76 


1.41 




sPP/PE 0.50 


0.0328 


13.89 


1.19 




iPP/PE 0.25 


0.0328 


11.35 


1.47 


0.005" 


iPP/PE 0.75 


0.0328 


11.33 


1.48 


0.01 " 


hhPP/PIB 0.50 


0.0343 


12.4 


1.28 


0.027- 11. l/r'-'''" 


^Reference 138). 
'■Reference I49|. 
'^Reference [13|. 

Reference [14]. 
•^Reference 1501. 



occupies a particular site is determined at random. Even 
though the molecular center-of-mass coordinates from a 
UA MD simulation can be used as an initial point, our 
calculations were started afresh as a more stringent test 
of our procedure and to allow us to increase the num- 
ber of particles in the system (or equivalently, the spa- 
tial dimension) at will, capturing relevant features of the 
effective pair interaction potential and to improve the 
statistical sampling of the ensemble. The potential and 
its corresponding derivative are entered as tabulated in- 
puts to the program, and linear interpolation was used 
to determine values not found in the supplied numeri- 
cal data as the algorithm proceeds. Each site is given 
an initial velocity pooled from a Mersenne Twister ran- 
dom number generator jST]. Subsequently, the system is 
evolved using a velocity Verlet integrator. Equilibrium is 
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FIG. 1: Comparison of the effective pair interaction potential Vai3{r) derived from the HNC closure for the hhPP/PE blend, 
</> = 0.5, with x/Xs G {0.0,0.1,0.3,0.5,0.7,0.9}. The upper panels show v'^'^{r) obtained via the Fade approximation and the 
lower panels show w'^'^(r) from the Debye form. The inset highlights the change in the repulsive part of the potential as the 
reduced temperature is changed. The solid line represents the athermal regime (x/Xa = 0.0). In both the AA and BB curves, 
the repulsive component decreases as the system approaches the spinodal {x/Xb = 1), whereas the AB curve increases. 



induced in the ensemble by rescaling the velocity at reg- 
ular intervals, and is manifested in the system when ob- 
serving a Maxwcll-Boltzmann distribution of velocities, a 
steady temperature, a stabilized Boltzmann _ff-theorem 
function, and a decayed translational order parameter. 

Once the equilibration step is completed, velocity 
rescaling is discontinued and trajectories are collected 
over a traversal of ^ 8Rg. The fastest benchmark in our 
studies is for a duration of ^ 4 hours for a system con- 
Extensive mesoscale simulations were performed on 
a typical system, the hhPP/PE mixture, to investi- 
gate the consistency of our approach. Simulations 
were performed for compositionally symmetric mix- 
tures, but also while approaching the spinodal, x = 
{0.008, 0.012, 0.016, 0.019}, while changing the fraction of 
A and B species in the melt such that (f> = {0.5, 0.7, 0.9}. 



sisting of ~ 6000 particles, performed on a single-CPU 
workstation. This compares extraordinarily well with a 
UA implementation that requires ~ 24 hours for a system 
with 1600 particles performed in parallel on a 64-node 
cluster |14j for an equivalent trajectory. We stress that 
our benchmarks for mesoscopic simulations represent an 
underestimate in the computational time since these have 
not been subjected to a parallelized algorithm. 

Mesoscale simulation parameters for all of the hhPP /PE 
systems are presented in Table [TTl For systems with x ap- 
proaching Xs, larger simulation boxes, with 10,648 par- 
ticles, were used to properly account for the increase in 
the lengthscale of concentration fluctuations. Those sys- 
tems also required longer equilibration. These simula- 
tions were run using the LONI TeraGrid system[52] to 
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TABLE II: Mesoscale Simulation Parameters for Blends of hhPP/PE 



Form Factor 


Interaction Parameter 


Particles 





L/2 [R- 


Pade 


X/X« = {0.1,0.3,0.5,0.7} 


5324 


0.5 


8.549 


Debye 


x/x. = {0.1,0.3,0.5} 


5324 


0.5 


8.549 


Debye 


X/Xs = 0.7 


10,648 


0.5 


10.771 



Debye x = {0.008, 0.012, 0.016, 0.019} 10,648 {0.5,0.7,0.9} 10.771 



facilitate performing numerous simulations at a time. 



V. TOTAL PAIR CORRELATION FUNCTIONS 
OF THE POLYMER MIXTURE FROM 
MESOSCALE SIMULATIONS 

From the trajectories of our mesoscopic simulations, 
the intermolecular total pcfs are computed. Initially, we 
set X = to determine the liquid structure far from the 
spinodal temperature, i.e. under athermal conditions, 
(1 — x/Xs) — 1. Mesoscale simulation parameters for 



these blends are presented in Table III For these simu- 
lations we compare the resulting pcfs to UA MD simu- 
lations. To obtain center of mass distribution functions 
from UA simulations, the center of mass coordinates for 
each chain are evaluated at each time step as the aver- 
aged sum of the position coordinates of each united atom. 
The radial distribution of center of mass coordinates is 
evaluated in the usual method employed for liquids. [55] 
The resulting pcfs are shown in Figure [2] for the systems 
listed in Table |l] Mesoscopic simulations are found to 
yield a coarse-grained liquid structure in agreement with 
our theoretical predictions from the analytical expression 
of Equation [9] serving as a self-consistent check of our 
determination of the effective pair potential through the 
HNC closure. The results presented in Figure [2] were ob- 

TABLE III: Mesoscale Simulation Parameters for Athermal 
Blends 



System 


Particles 


<!> 


L/2 [R-^] 


hhPP/PE 


5324 


0.5 


8.549 


PIB/PE 


4096 


0.5 


8.365 


PIB/sPP 


5488 


0.5 


10.416 


sPP/PE 


1728 


0.5 


5.635 


iPP/PE 


4913 


0.25 


8.482 


iPP/PE 


1728 


0.75 


6.016 



tained using the Fade approximation (Equation |6| which 
works sufficiently well under athermal conditions where 
the low k behavior is less important since critical fluctu- 
ations are assumed to be small. 

The liquid structure from mesoscopic simulations are 
in general consistent with data obtained from UA MD 
simulations, with the exception of blends containing iPP 
and PIB for which theoretical predictions and meso- 
scopic MD predict a less pronounced correlation hole 



than UA MD simulations. These observations are not 
surprising since these systems tend to possess very efh- 
cient intramolecular packing, leading to smaller isother- 
mal compressibilities and thermal expansion coefficients 
when compared to other polyolefin blends [TU [SD] . The 
effective intramolecular packing arises from the attractive 
interactions between methyl moieties induced by their ge- 
ometrical arrangement. However, the theory and meso- 
scopic simulations do exhibit good agreement for r Rg. 

Moving to the thermal regime, where large scale fluc- 
tuations in the local concentration develop as the system 
approaches a second order phase transition, we present 
results for the typical 50:50 mixture of hhPP/PE, al- 
though the theory and methods employed are ubiquitous 
and generally applicable to a wide range of systems. For 
these simulations the value of the x parameter was varied 
such that x/Xs = {0.1,0.3,0.5,0.7}, in order to see the 
changes in the pcfs as the system approaches the spin- 
odal. Figure |3] shows the dependence of the partial corre- 
lation functions on the interaction parameter, x- The left 
side of Figure [3] shows the resulting correlation function 
from mesoscale simulations with the potential obtained 
via the Fade approximation, after using our truncation 
scheme for c'^'^(r) in the HNC (see upper panels Fig [I]). 
Use of the Fade approximation has the advantage of al- 
lowing a fully analytic solution for h'^'^{r), as shown in 
Equation [9] which shows quantitative agreement with 
UA simulations in the athcrmal limit (as shown in Figure 

The right panels of Figure |3] shows the correlation func- 
tion obtained using the potential derived from the Debye 
form (see lower panels of Fig[l]). Here, comparison is 
again made to numerical predictions based on Equation 
3 since an analytic solution, such as that of Equation 
9 is not possible when the Debye form is used. In both 
the right and left panels of Figure |3] mesoscopic simu- 
lations show quantitative agreement with our theoreti- 
cal predictions, indicating the self-consistency of our ap- 
proach. Furthermore, despite the differences in the po- 
tential used in the simulation, Figure |3] shows that the 
resulting pcfs from either the Fade or Debye form are 
qualitatively similar. Lastly, we note that despite the 
approximations made in obtaining the analytical form of 
Equation |9j our analytical expressions recover the cor- 
rect k=0 limit. [331134] In fact, all of the forms for h'=''{k) 
exhibit the same k = behavior. 

The standard approach to describe the mixing behav- 
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FIG. 2: Comparison of mesoscopic simulations [open symbols] with UA MD simulations [filled symbols] for the hai3{r) of polymer 
mixtures under athermal conditions. Also shown are theoretical predictions [solid curves] based on our analytic expression, 
Equation[9] Presented are data from AA [circles], AB [triangles], and BB [squares] contributions for compositionally symmetric 
and asymmetric systems. For comparison, numerical predictions obtained from Equation [S] using the Debye form are shown 
[dashed curves]. For clarity the inset highlights the peak region for each separate contribution. 



ior of polymers is the Flory-Huggins model. Under Flory- 
Huggins treatment, the phenomenon of demixing is un- 
derstood in terms of contributions to the free energy of 
mixing. Generally, at low enough temperatures the trans- 
lational entropy, which is associated with the center of 
mass motion of the molecules and always favors mixing. 



is outweighed by local monomer-monomer interactions. 
In most cases, van-der-Waals interactions are stronger 
between like pairs than those between unlike pairs, re- 
sulting in a positive free energy of mixing. As a result, 
lower temperature favors spontaneous demixing due to 
changes in the local free energy of the system. [T] In an 
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FIG. 3: Comparison of mesoscopic simulations [symbols] with numerical predictions [curves] for the ha/3{r) of a 50:50 mixture 
of hhPP/PE for different values of the ratio, x/Xs- The left panel shows results obtained using the Fade approximation with 
our truncation scheme. The right panel depicts the results when the Debye form is used. Mesoscale simulations are shown to 
capture the structural changes that occur as the system approaches the spinodal. The inset highlights the peak region of h{r) 



empirical manner, the Florry-Huggins parameter, Xi is 
used to describe these changes in local free energy. At 
the hmit of the spinodal temperature, X ~^ Xsj smd since 
X ot: ^, positive values of x always lead to incompatibility 
of the mixture. [T] 

In real systems, the simple Flory-Huggins model does 
not hold, and the x parameter may be a complicated 
function of N, </>, and T, leading to the variety of phase 
behaviors observed in polymer blends. For example, 
some blends phase separate upon cooling, while others 



show an opposite trend in demixing and phase separate 
upon heating. It is customary to fit the experimen- 
tal temperature dependence of a mixture to the form 
X — a + b/T where a and b may be either positive or 
negative depending on the system. Table [l] shows the 
experimentally determined a and b parameters for a few 
of the systems investigated in this paper. It should be 
noted that when applying an equation for the x param- 
eter from the literature, the x value must be normalized 
by the average number of UA sites per monomer [14] to be 
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consistent with the site-basis description adopted here. 

In our present treatment, the interaction parameter, x 
is treated as an adjustable parameter, which describes the 
interactions that drive phase separation. It is analogous 
to the Flory-Huggins parameter; however, since in our 
model it represents a system specific parameter, it may 
be given any value positive or negative depending on the 
behavior of the system of interest. The advantage of 
a mesoscale approach is that once the system specific 
parameters are defined, the trends in phase behavior can 
be readily calculated without requiring restrictively large 
MD simulations. 

As a further implementation of our theory, we per- 
form mesoscale simulations at several fixed values of x 
for which the fraction of A and B species in the melt is 
varied. For these simulations, we again use hhPP/PE as 
a typical system and vary the volume fraction such that 
(p = {0.5,0.7,0.9}. In order to better capture the large 
scale structural changes, simulations were performed in 
a large box with 10,648 particles. Figure [4] shows the 
resulting pair correlation functions for mesoscale simula- 
tions run with x = 0.008 and x = 0.012, and Figure [5] 
shows the case where x = 0.016 and x = 0.019. In all 
cases, mesoscale simulations correctly capture the struc- 
tural changes that depend on the concentrations of the 
species in the mixture when comparison is made with our 
theoretical predictions. For these simulations we limit 
our consideration to using only the Debye form in Equa- 
tion[3]to avoid any effects due to the truncation scheme in 
the low k region. Once more, theory and mesoscale sim- 
ulations appear to be fully consistent in predicting the 
structural information of the mixture in the lengthscales 
larger or equal to the polymer radius-of-gyration. 



VI. SCATTERING FUNCTIONS AND 
CONCENTRATION FLUCTUATIONS 

The mesoscale pair correlation functions effectively de- 
scribe the polymer fluid as a liquid of soft colloidal par- 
ticles. Once these pcfs are obtained from simulation, 
any property of the liquid can be calculated, including 
the equation of state, internal energy, compressibility, 
and others. [45] In this section, we examine the extent to 
which our classical MD simulations of soft colloidal par- 
ticles reproduce the structural changes which occur as 
the system approaches the spinodal. Due to the increas- 
ing length scale of fluctuations as the system approaches 
the critical temperature for demixing, UA simulations 
can only reach a very limited region of the phase dia- 
gram. An advantage of using a procedure that captures 
the structure at the mesoscopic scale is that the relevant 
length scale of the simulation can increase considerably 
with respect to UA MD, and simulations can describe the 
increasing lengthscale of the fluctuations. Thus a meso- 
scopic picture greatly facilitates the ability to capture 
this phenomenon, since we are able to simulate many 
thousands of chains represented as soft spheres. Mod- 



els using Monte Carlo methods with phenomenological 
potentials have been previously performed at the level 
of soft colloids, demonstrating the valuable information 
that may be gained about phase transitions. [31] [SI] The 
advantage of the procedure presented here is that the 
potentials used to simulate the system are explicitly pa- 
rameter dependent, being related to the system-specific 
molecular parameters, such as Rg. The potentials ob- 
tained in this manner allow for mesoscale simulations to 
be performed on any number of different, but specific, 
systems under different thermodynamic conditions, map- 
ping them as soft colloids. 

The static structure factors for each component are 
calculated from our simulations by Fourier transform of 
the total correlation function, 

/■°° sin kr 

SAAik) = (j) + ■iTTcjP pch / r"^ —^hAA{r)dr , 
f°° sinfcr 

SBB{k) = 1 - 47r(l - (jyfpch / r^ ——hBB{r)dr, 

Jo 

SAB(k) = 4^,^(1 - cb)p,h I r^^-^hAB{r)dr{16) 

Density and concentration fluctuation contributions can 
be written as linear combinations of the static struc- 
ture factors according to the formalism of Bhatia and 
Thornton|36). Here, the density fluctuation, S'''''(fc) is 
given by 

SPP{k) = SAAik) + SBBik) + 2SAB{k) . (17) 

The concentration fluctuation contribution, S'^'^{k) may 
be expressed as 

S^^ik) = {l-^fSAA{k) + cl)^SBB{k)~2cj){l~cl))SAB{k) , 

(18) 

and is particularly important since it provides informa- 
tion about the stability of the binary mixture against 
demixing. The coupling term, SP'^{k), is given by 

SP^'ik) = (1 - cf,)SAA{k) - <l>SBB{k) + (1 - 2cf,)SAB{k). 

(19) 

Figure [6] shows the colloidal partial structure factors, 
SPP(k), S't"t'(k), SP't'ik), calculated from pcfs obtained 
from mesoscopic simulations shown in the right panel of 
Figure [3] using Equations [T6][l9] The data from the simu- 
lation is compared to predictions based on our numerical 
values for /i'"^(fc), obtained from Equation [s] using the 
Debye function. Since it is particularly pertinent to cap- 
ture the low k behavior where concentration fluctuations 
will diverge as the spinodal is approached, we use the 
results for the Debye form since the Fade approximation 
introduces unphysical effects in this regime, typically for 
kRg < 2. As seen in Figure [6| the curves of the den- 
sity fluctuation contribution, S^P{k), which behaves sim- 
ilarly to the static structure factor for a single-component 
liquid. [34] are indistinguishable over the range of x 
vestigated. The function S'"^{k) exhibits a slight depen- 
dence on the ratio x/Xs in which the minimum at low 
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FIG. 4: Comparison of mesoscopic simulations [symbols] with numerical predictions [curves] for the hap{r) of hhPP/PE for 
different values of 0. Left panels show data when x = 0.008. Right panels show data for x = 0.012. Shown are the separate 
contributions for AA [circles], AB [triangles], and BB [squares] interactions. As <j> increases, the fraction of species B in the 
simulation box decreases, and thus, the statistics become poorer for BB interactions. 



k becomes slightly more pronounced. The minimum in 
S'"^{k) represents the length scale for asymmetry in the 
mixture arising from the difference in particle size.|34j 
The partial structure factor, S'^'^{k), exhibits a charac- 
teristic diverging behavior as the spinodal is approached, 
indicating an increase in the length scale of concentration 
fluctuations. 

As illustrated in the upper left of Figure |6) S"^'^(0) 
increases as the ratio x/Xs —^1. As the system nears 



the phase transition, the divergence of S'^'^{k) is indica- 
tive of the concentration fluctuations becoming increas- 
ingly macroscopic. Since concentration fluctuations oc- 
cur on an increasingly large scale, the relevant region of 
the S'^'^{k) curve occurs in the low-k region; hovifever, 
due to periodic boundary conditions, simulation data is 
only reliable at a distance less than half the length of the 
simulation box. This makes extrapolation of the k=0 
limit from mesoscopic simulations still difficult, as seen 
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FIG. 5: Same as Figure [4] except that left panels show data when x = 0.016, and right panels show data for x ~ 0.019. 



in Figure [6j even though thousands of particles were rep- 
resented. In this respect, our numerical predictions may 
serve as a guide for extending S{k) to the k=0 limit. 
Furthermore, we have previously shown that Equation [9j 
also gives an estimate for S"^'^(0) given bv[5I] 



1- 



^ii-m.^-ire, ^^^^ 



Even though it is based on the Fade approximation. 
Equation 20 may be used to estimate ^"^"^(O) since h'^'^{k) 



calculated from the Fade approximation has the same 



k — Q limit as h'^^{k) from the Debye form. The lower 
right of Figure [g] shows a linear plot of 1/ S'^'^iO) vs. x for 
which the fc = limit was determined by our theoretical 
predictions. 



Following Equations [T6][T9) the concentration fluctua- 
tion partial structure factor, S'^'^{k), was calculated from 
the mesoscale simulations presented in Figures |4] and [5] 
where the volume fraction, 0, was changed. The result- 
ing S"^'^(fc) is presented in Figure [t] along with theoret- 
ically predicted values using the Debye formula. Once 
again, mesoscale simulations show an increase in concen- 
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FIG. 6: Top Left: Partial structure factor, 5"*'^(fc), obtained from mesoscopic simulations [symbols] of the coarse-grained 
mixture of 50:50 hhPP/PE with x/Xs & {0.0,0.1,0.3,0.5,0.7}. The curves represent theoretical values obtained using the 
Debye function. Top Right: Partial structure factor, S''''{k) and Bottom Left: S'"^{k) are also shown for different values 
of x/Xs- S''''{k) does not change noticeably with x/Xs but S'"^{k) has a slight x/Xs dependence at low k. Bottom Right: 
Extrapolated 1/S"^'*(0) values vs. x [symbols]. The line represents a linear fit to the data and is extrapolated to the spinodal, 
Xs (dashed line). 



tration fluctuations as the thermodynamic conditions are 
changed, and x ^ Xs or (f) 0.5. In general, mesoscale 
simulations are consistent with our theoretical predic- 
tions based on Equation [3] up to the limit set by the 
finite box size. As seen in Figure [7| when x is low or 
the polymer volume fraction of one species is large, the 
system is well mixed and the extrapolation to low k is 
straightforward. However, for the case when (j) = 0.5 and 
X = 0.019, as depicted in the lower right panel of Fig- 
ure [Tj it becomes more difficult to reach the fc = limit 
from mesoscale simulation, even if the precision is higher 
than for atomistic simulations for the reasons previously 
discussed. Since our simulations are consistent with our 
theoretical predictions as shown in Figures [2] - [5j we es- 
timate the extrapolated S{k = 0) limit based on these 
predictions. 

Once this method is employed, it is possible to discern 
the phase behavior of the mixture from the extrapolated 
k ~ limit. In order to include more data points, we 
calculate ^"^"^(O) for a range of x and (j) values, based 
on our solution to Equation |3] These are presented in 



Figure [5] which shows the structure factor as a function 
of the volume fraction for several fixed values of the x 
parameter. The interpolation between the points is given 
by Equation |20[ which demonstrates that our analytical 
expression is useful in determining the phase behavior. 

Finally, in Figure [9] a plot of the inverse structure fac- 
tor, S''^'^(0), vs. X each value of (j) shows the lin- 
ear behavior from which the spinodal, Xsi may be ex- 
trapolated and used to sketch the phase diagram of the 
system. In the bottom panel of Figure [9] the spinodal 
curve is compared to the predicted Flory-Huggins model, 
Xs = [2Af^(/)]~-^-|-[2A'^B(l—0)]~^, which was used in Equa- 
tion[8j The spinodal curve from our simulation exhibits a 
characteristic parabolic shape consistent with mean-field 
theory, where ''^ — xlXs)~'^ i v = 1/2. In the imme- 
diate region of the critical temperature, mean-field theory 
breaks down, and Ising-type critical behavior is expected. 
For this narrow temperature region, the linear extrapo- 
lation in Figure |9] would be invalid and the spinodal will 
exhibit a flatter peak. [M] For long polymer chains, the 
temperature region for which mean field theory becomes 
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FIG. 7: The concentration fluctuation partial structure factor, 5"'"*(fc), calculated from mesoscale simulations [filled symbols] 
at different values of <j) for the mixture hhPP/PE. The curves represent theoretical predictions. 




FIG. 8; The extrapolated it = limit of S"*''(fc) based on 
our numerical predictions [data points] and from our analyt- 
ical expression, Equation 20 [curves] as a function of for 



different fixed values of x, for the mixture hhPP/PE. 



invalid is very small, since the Ginzburg number, which 
scales inversely with N, is small. [TT] As seen in the upper 
panel of Figure [9] most of the simulations performed are 
well within the temperature region described by mean- 
field theory. Although the linear extrapolation becomes 
less quantitative near the horizontal axis, the mean-field 
approximation is consistent with our data. 



VII. CORRECTIONS TO THE DEBYE 
INTRAMOLECULAR FORM FACTOR 

Upon examination of Figure [2] it appears that there 
is slightly better agreement between UA MD simula- 
tions when compared to our analytical results using the 
Fade approximation than with the full Debye form (as 
indicated by the dashed lines). Since the Fade form is 
an approximation, this improvement is likely due to a 
cancellation of errors. The Debye formula is exact for 
ideal Gaussian chains; however, Wittmer, et. al, have 
recently shown that dense polymer melts exhibit devi- 
ations from ideal Gaussian behavior because of long- 
range correlations arising from the repulsive interaction 
of chain segments. (55j These deviations become more sig- 



16 




= 0.5 
= 0.6 
= 0.7 
= 0.75 
= 0.8 
= 0.85 



[2TI was input into Equation [3] to obtain a corrected form 
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FIG. 9: Top: Inverse concentration fluctuation structure fac- 
tor, S"*'^(0) plotted against the interaction parameter, x foi" 
diflerent values of </!>, for the mixture hhPP/PE. The solid 
line depicts a linear fit to the data and the dashed line shows 
the extrapolation to the spinodal. Bottom: Phase diagram 
for the coarse-grained mixture obtained from the above ex- 
trapolation to the spinodal, Xs- The solid curve depicts the 
Florry-Huggins analytical expression. 



nificant for polymers confined between walls in ultrathin 
films. [56] In this section we investigate the implemen- 
tation of higher order corrections to the Gaussian ap- 
proximation on the effective pair potential by evaluating 
Equation [3] numerically with a corrected from of the in- 
tramolecular form factor. 

In the infinite chain limit (N — >■ 00) it has been pro- 
posed that corrections to the Debye formula in the inter- 
mediate wave vector range depends only on the monomer 
density, such that [551 [57] 



1 



1 



1 fc3 



a;™™(/c) u;^,!^^,(fc) 32 p 



(21) 



Although approximate for finite chain lengths. Equation 



of the pair potential which is shown in Figure 10 (top 
left) for a 50:50 mixture of hhPP/PE. The resulting cor- 
relation functions, displayed in Figure [TO] show that the 
corrected Debye formula agrees very well with UA MD 
simulations for this sample, indicating that the disagree- 
ment between mesoscale simulations using the Debye for- 
mula and UA simulations on intermediate lengthscales is 
due to non-Gaussian behavior of real chains as the Flory 
ideality hypothesis breaks down. On the local scale, how- 
ever, the corrected Debye and the UA-MD simulations 
tend to disagree. This is not relevant for systems with 
long chains, such as the hhPP:PE mixture, but it be- 
comes important for short chains, e.g. mixtures of PIB 
chains, where the behavior at short distance becomes un- 
physical. In conclusion, while in the current publication 
we limit our investigation to just this correction term for 
the hhPP:PE mixture, further study is necessary to inves- 
tigate if the observed improvement is a common feature 
of long-chain mixtures, independent of their monomeric 
structures. The pcfs obtained using the Pade approxima- 
tion (Figure [2]) are also shown in Figure 10 and compare 
well with the corrected Debye results. 



VIII. APPLICATIONS TO MISCIBLE LCST 
BLENDS 

While most polymer blends are immiscible and tend 
to demix at experimentally relevant temperatures, some 
systems are known to be miscible having a lower critical 
solution temperature (LCST). In this section we demon- 
strate the extension of our approach to model LCST 
blends where the effective x parameter may be nega- 
tive over most of the temperature range of interest. It 
is worth noticing that while the hhPP/PIB blend is mis- 
cible, the iPP/PIB blend is immiscible, indicating that 
subtle changes in the specific polyolefin architecture may 
give rise to a completely different phase diagram. The 
temperature dependence of the x parameter for the mis- 
cible hhPP/PIB blend is reported in Table |l) The x pa- 
rameter in the literature is defined on a monomer ba- 
sis and must be divided by the number of united atom 
sites per monomer (4.8 for hhPP/PIB) to be consistent 
with the UA site description used here. We performed 
mesoscale simulations for various temperatures of a mix- 
ture of 50:50 hhPP/PIB {xs = 0.021). 

The resulting correlation functions determined for two 
temperatures, 2000K and 200K, from mesoscale simula- 
tions are shown in Figure [TT] When compared with Fig- 
ure [3] it is evident that the pcfs for the hhPP/PIB blend 
exhibit an opposite trend with temperature. These differ- 
ences are clearly evident in the concentration fluctuation 
structure factor, which was calculated from these pcfs at 
various temperatures and is shown in the bottom right of 
Figure[TTJ As depicted in the low wave vector behavior of 
S'^'^{k), fluctuations in the concentration become smaller 
as the temperature is decreased, and the system becomes 
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FIG. 10: Top Left: The effective pair potential between A type chains for an athermal mixture of hhPP/PE with = 0.5 
when corrections to the Debye formula are included (solid line). The dashed line indicates the potential obtained using the 

uncorrected Debye forrmila. Top Right: The A A component of the correlation function calculated from mesoscale simulations 
using the corrected Debye fornmla at x — 0.0 (open circles) and x ~ 0.7 (open triangles). The solid line represents theoretical 
predictions, and the dashed line indicates predictions using the Debye formula. Filled circles represent UA MD simulations. 
Bottom Left: BB component and Bottom Right: AB component of the correlation function for the same mixture. The pcfs 
obtained using the Fade approximation are shown to nearly superimpose with the corrected Debye form. (Partially shaded 
circles). 



more stable. These results indicate that our procedure of 
mapping polymer blends as soft-colloids and performing 
nicsoscopic simiilations using an effective pair potential 
can be applied to miscible LCST blends given that the 
temperature dependence of the x parameter is known. 



IX. SUMMARY 

In this paper, we present the implementation of our an- 
alytical coarse-grained description for polymer mixtures 
in mcsoscopic modeling through computer simulation. 
Using the analytical representations for the intermolecu- 
lar total pair correlation functions at the center-of-mass 
level and the hypcrnettcd-chain closure, we derive the ef- 
fective pair interaction potentials which are the required 



input to the coarse-grained simulations. The simulations 
are then carried out and the coarse-grained liquid struc- 
ture, as given by the intcrmolecular pair correlation func- 
tion, is extracted from the trajectories. In the athermal 
regime, results are compared with our theoretical pre- 
dictions and data obtained from united atom molecular 
dynamics simulations. In the thermal regime, mesoscopic 
simulations capture the relevant trends for demixing of 
polymers in the miscible regime approaching the spin- 
odal. These results are used to calculate static structure 
factors which are related to the increasing concentration 
fluctuations of the mixture. By extrapolation to the low 
wave vector limit, we are able to determine the phase di- 
agram of the coarse-grained mixture which is consistent 
with mean-field theory predictions. The consistency of 
all representations supports the theoretical foundation of 
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FIG. 11: Top Left: AA component of the correlation function for the miscible blend, hhPP/PIB (0 = 0.5) at T = 2000iC 
(circles) and T = 2QQK (triagles). Theoretical predictions are indicated as solid lines. Top Right: BB component and Bottom 
Left: AB component of the same mixture. Bottom Right: The concentration fluctuation structure factor for hhPP/PIB 
obtained from mesoscale simulation (symbols) and from theory (solid line) at various temperatures. 



our development. 

While other methods exist to obtain equilibrium prop- 
erties of blends, an analytical potential is desirable for 
many of these approaches. For, example, the Gaussian 
equivalent representation (GER) has been implemented 
using a purely repulsive Gaussian potential for a sys- 
tem of interacting particles to obtain a field-theoretic 
representation of the partition function, and used to 
compute structural and thermodynamic quantities of 
interest. [581 ES] Possible future applications of the ana- 
lytical potential derived in this publication could include 
using it in low-cost approximation methods, such as the 
zeroth-order GER model formalism. 

Although the current publication is focused on deter- 
mining the equilibrium structure of blends under various 
conditions, the proposed procedure of obtaining the ef- 
fective potential and performing mesoscale simulations 
should be useful in determining the non-equilibrium and 



dynamic behavior of these systems as well, where time 
must be treated explicitly. Future work in this direc- 
tion should include examining the phase transition that 
occurs when the thermodynamic conditions of the sys- 
tem are suddenly changed. In the context of a multiscale 
modeling procedure, mesoscale simulations, such as those 
performed here, may be coupled to a more detailed de- 
scription in order to combine local and global information 
over multiple length and time scales. 
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